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I.  INTRODUCTION 

A forward  Monte  Carlo  method  has  been  developed  to  solve  the  time-dependent 
heat  conduction  or  diffusion  equation.  The  general  method  has  been  implemented 
so  as  to  cover  a variety  of  boundary  conditions. 

The  method  is  based  on  a "floating  volume"  random  walk,  similar  to  the  one 
in  eui  adjoint  method^.  A novel  problem  has  been  posed,  and  successfully  solved, 
on  the  treatment  of  boundary  conditions.  The  solution  requires  the  introduction 
of  a biasing  function,  leading  to  a biased  random  walk.  The  first  step  of  the 
random  walk  is  particular  to  the  forward  method.  The  succeeding  steps  can  be 
considered  self-adjoint,  as  they  are  identical  in  the  forward  and  adjoint  case. 

It  is  hoped  that  our  study  of  biased  random  walks  will  also  prove  useful  in  any 
future  development  of  importcince  biased  adjoint  methods. 

The  detailed  analysis  and  computer  implementation  are  still  under  way.  The 
details  of  sampling  algorithms  have  been  completely  worked  out  in  the  case  of  known 
temperature  boundary  conditions.  Sections  VIII. 2 through  IX. 2. 2,  as  well  as 
Appendix  B,  are  confined  to  this  type  of  boundary  condition.  The  remainder  of  the 
report,  including  Appendix  A,  applies  to  general  "radiation  type"  boundary  con- 
ditions. 
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II.  DERIVATION  OF  AN  INTEGRAL  EQUATION 


Let  us  consider  the  heat  conduction  equation 


DV  ^ T(x,t)  - - 0 

X d C 


(1) 


defined  over  a volume  enclosed  on  a surface  with  the  initial  conditions 


T(x,0)  given  for  xefl. 


(2) 


and  the  boundary  condition 


a(x)  ° Tg(x,t)  - T(x,t) 


(3) 


T (x,t)  given 

g 


The  problem  is  to  estimate  the  temperature  profile  at  a given  time  t^.  For 
that  purpose,  it  will  be  useful  to  introduce  three  Green's  functions  G^,  i=0, 

1,2  satisfying  the  following  equation: 


D G^(x,x' ,tQ-t)  + G^(x,x',tQ-t)  * 0,  xea 


(4) 


with  the  initial  condition 


G^(x,x',0)  » 6(x-x')  x',xen^ 
and  the  boundary  condition 

G^{x,x',tQ-t)  » - a(x)  Gj|^(x,x',tQ-t)  x',xe4l 
a(x)_>P  ^^—^0 


(5) 


(6) 


* 

For  simplicity  of  discussion,  we  assume  that  the  diffusion  coefficient  l>K/pc 
is  constant.  The  method  is  readily  generalized  to  the  case  of  D^constant  in 
finite  geometrical  regions,  with  discontinuous  variation  of  K,  p,  and  c across 
boiindaries. 
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Multiplying  Equation  (1)  by  (x,x' jt^-t) , Equation  (4),  written  for  i*l,  by 
T(x,t),  subtracting  and  integrating  over  and  over  time,  we  obtain,  after 
applying  Green's  theorem 


T(x,t)  D Gi(x,x',tQ-t)J  dS^ 
+ T(x,t)  Gj^(x,x',tQ-t)”j  dV^ 


Taking  into  account  the  boundary  conditions  (3)  and  (6)  is  the  surface  term, 
and  the  initial  conditions  (2)  and  (5)  in  the  volume  term,  we  obtain: 


D ^Gi(x,x',to-t)dSx 


G (x,x',t-)T(x,0)dV„  - T(x',t„) 


0 


or,  interchanging  x and  x*: 


T(x,tQ)  = V(x)  + S(x) 


(7) 


where 


S(x) 


T(x',0)  G^  (x',x,tQ 

'fyk 


° Tg(x',tQ-t)  ^G^{X',X,t)dS^, 


(8) 


(9) 


The  object  of  forward  Monte  Carlo  is  to  generate  a population  of  weighted  points 
xefjj^  with  a density  for  given  t^. 
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A simple  and  correct  Honte  Carlo  approach  to  this  problem  is  to  consider 


Equation  (7)  as  the  sum  of  two  terms  to  be  sampled  separately. 

The  volume  terms  Equation  (8)  presents  ikj  difficulty.  One  can  Scunple  x* 
from  a prob^dJility  distribution  function  (pdf)  proportional  to  T(x*,0).  Given 
x'  and  t^,  one  can  sample  x from  (x',x,tQ)  using  the  self-adjoint  method  de- 
veloped previously. 

The  surface  term  Equation  (9)  is  not  as  straight  forw^u:d  to  handle.  The 
kernel  - — p G {x',x,t)  cannot  be  reduced  to  a pdf  because 


f 


^ * t-vO 


Let  us,  however,  define  a function  G2(x,x',t)  which  satisfies  Equations  (4-6) 
for  i=2,  and  the  function 


Q(x',t) 


G2(x',x,t) 


dV 


X 


(10) 


It  is  shown  in  Appendix  A that  if  the  surfaces  and  has  the  same  outward 
normal  at  x'  (see  Figure  1) , then 


which  implies  that  the  kernel 


T (x’,t  -t)Q(x',t) 
g 0 

and,  given  x'  and  t,  sample  x from  a pdf  proportional  to  the  kernel  (11). 

The  function  Q(x',t)  can  be  considered  as  a biasing,  or  importance  function. 
It  is  defined  (Equations  10  and  4)  over  a "floating  importance  volume" 
which  we  define  more  precisely  as  the  superscript  x'  expressing  the  fact 
that  the  surface  is  tangent  to  at  x' . The  shape  of  is  otherwise 
arbitrary;  it  can  be  chosen  simply  enough  such  that  Q be  analytically  known.  In 

X * 

practice,  we  will  eventually  restrict  the  surface  to  be  an  infinite  plane 
tangent  to  2)^  at  x'. 

As  long  as  we  have  the  necessity  to  introduce  a biasing  function  for  samp- 
ling the  surface  term  (9) , we  find  it  also  desirable  to  introduce  a biasing 
function  for  the  volume  term  (8),  In  order  to  define  such  a function,  let  us 
first  generalize  the  definition  of  "floating  importance  volume".  The  floating 
volume  ^2  been  defined  for  points  x'  on  Let  us  now  consider  a point  y 

internal  to  the  configuration.  Let  x'  be  the  point  on  closest  to  y (see 
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Figure  1).  we  define  the  "floating  importance  volume"  associated  with  y as 
any  volume  its  surface  being  tangent  to  at  x'.  We  generalize  our 
definition  of  superscripts  and  denote  that  volume  as  The  restricted  de- 

finition of  the  last  paragraph  is  preserved  as  y and  x'  coincide  if  y approaches 
the  surface  2^.  Let  us  now  define  the  function 


Ey(z,t)  = 


(z,x,t)dV^ 


_ y 

where  G^(z,x,t)  satisfies  Equations  (4-6)  for  i=2  and 

We  propose  to  useEq.Ey (x,t)  as  a biasing  function  for  internal  points 
y,  and  rewrite  the  volume  term  (8)  in  the  form 


V(x)  = 


G (x' »x,t  ) 

^ ..  dv 


IV.  A PRACTICAL  MONTE  CARIO  ALGORITHM 


Equation  (14)  can  be  rewritten  as 


V(x)  = T^Uy{x) 


(15) 


where 


r 


\ T(x*,0)E^(x',tQ)dV^, 


(16) 


and 


0,  ,,  Gi<’‘''X,tQ) 


E;;:(x--,-t")- 


(17) 


p^°(x')  = T(x*,0)E^. (x',tQ)/T^ 


(18) 


Similarly,  Equation  (12)  can  be  rewritten  as 


S(x)  = TgUg(x) 


(19) 


where 


Tg  = / dtjj-  DTg(x*,tQ-t)Q(x',t)  dS^, 


(20) 


0,.  / 0,  . . . " Gj^(x',x,t) 


L Ps  Q(x-,t) 


dS  , (21) 


Pg°(x',t)  = DTg(x',tQ-t)Q(x’,t)/Tg 


(22) 


Substituting  (15)  and  (19)  into  (7) , we  obtain 


T(x,tQ)  = T^U^(x) 


(23) 


where 


T-  “ T,.  + T„ 
TVS 


(24) 
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and 


U^(x)  - PyU^(x)  + {l-P^)Ug(x)  (25) 

Pv  ■ V^T 

As  stated  previously,  the  object  of  a foi-ward  Monte  Carlo  method  is  to  gen- 
erate a population  of  weighted  points  representing  the  density  TCxjt^)  for  given 
tp.  We  define  a random  member  of  such  a population  as  a sample  of  TCxjt^),  and 
the  process  of  generation  of  such  a member  as  seunpling  TCXft^). 

According  to  Equation  (23),  in  order  to  sample  T(x,tQ),  we  can  sample  U^(x) 
and  multiply  the  weight  by  To  sample  U^(x),  consider  Equation  (25).  With 

probeibility  p,^,  sample  U^(x),  else  sample  Ug(x)„ 

We  now  turn  to  prescribe  methods  to  san^le  U^(x)  emd  Ug(x).  In  the  case  of 
U (x),  consider  its  definition  (Equation  17).  To  szunple,  we  can  first  sample  x' 
from  the  pdf  p^^(x')  defined  by  Equation  (18)  and,  given  x’  (£md  t^) , s^unple  x 
from  the  kernel 

K^(x*,x,tQ)  = Gj^(x',x,tQ)/E^(x',tQ)  (27) 

In  the  case  of  Ug(x),  consider  its  definition  (Equation  21).  To  saiT5)le,  we 
' can  first  san^le  x'  and  t from  the  pdf  Pg^(x',t)  defined  by  Equation  (22)  and, 

given  x*  and  t,  sample  x from  the  kernel 

' Kg(x',x,t)  - Gj^(x',x,t)/Q(x',t)  (20) 

The  sampling  of  the  kernels  (27)  and  (28)  can  be  achieved  by  constructing 
a random  walk  which  we  are  about  to  describe.  The  walks  for  and  are  identi- 
cal except  for  the  first  step. 


I 

i 

% 


THE  BIASED  FLOATING  VOLUME  RANDOM  WALK 


► 
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V. 1 The  Integral  Equation 

As  defined  previously,  the  configuration  under  investigation  is  bounded  by 
a surface  with  a volume  We  also  defined  a "sin^jle"  surface  with  a 

volume  Sl^,  tangent  to  at  a variable  point  x*  (see  Figure  2).  In  addition, 
let  us  define  another  "simple"  surface  with  a volume  which  is  wholly 

K 

contained  in  both  and  and  a Green's  function  G^  defined  over  which 
satisfies  Equations  (4)  eind  (6)  for  i«0.  We  rewrite  these  equations  euid  write 
the  initial  conditions  in  the  form 


Gq(x,x" 

,t) 

- |^Go(x,x",t)  . 0 

xcOo 

(29) 

Gq(x,x" 

,0) 

= 6(x-x") 

xefio 

(30) 

Gq(x,x" 

,t) 

= - 0i(x)  — GQ(x,x",t) 

xtfio 

(31) 

Let  ^0  “ ®1  ^ common  part, 

is  the  rest  of  which  is  Internal  to  and 
Let 

Oq(x)  » Oj^Cx)  - a^{x)  for  xeS^ 

and 

“ 0 xeSj^ 


FIGURE  3. 


if  any,  of  and 


(32) 

(33) 
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Multiplying  Equation  (4)  by  G^,  Equation  (29)  by  subtracting,  and  in- 
tegrating over  we  obtain,  after  applying  Green's  theoremt 

£ Gj(x,xi.tj-t)  - Gj(x,x',t^-t)  Ij- G(|(x,x-,t)]  dSj, 

[=o‘x-x-.t) 

+ G^(x,x',tp-t)  GQ(x,x",t)J  (34) 

The  Zq  “ integral  can  be  brolcen  up  into  an  integral  and  an  S^-integral. 
The  S^-integral  vanishes  because  of  the  boundary  conditions  (6),  (31),  and  (32). 
The  first  term  of  the  Sj^-integral  vanishes  because  of  the  boundary  conditions 
(31)  ^tnd  (33).  The  integrand  of  the  volume  term  can  be  rearrwged.  Equation  (34) 
can  be  rewritten  as 


/ ^GQ(x,X-,t)dS^ 

.|,/„^G,x.x..V.G„,x.x.. 


t)dV  - 0 
X 


Integrating  over  time  from  0 to  t^,  and  taJclng  the  initial  conditions  (5) 


and  (30)  into  account,  we  obtain: 


G^(x",x',tQ)  - GQ(x',x",tQ) 


10 


or,  dividing  by  E^„(x",tQ)  defined  by  Equation  (13): 


G^(x",x',tQ) 


Go(x',x”,to) 


-'0  s. 


G^(x,x',tQ-t) 

Ejj(x,tQ-t) 

Ex..(x,tQ-t)  ^GQ(x,x"t) 

E,(».to-t) 

- — 

(37) 


dS 


Equation  (37),  written  for  i»l,  can  be  considered  as  an  integral  equation  for 
K^(x",x',tQ)  defined  by  Equation  (27). 

Writing  Equation  (37)  for  i*2,  integrating  x'  over  and  talcing  (13)  into 
account,  we  obtain 


f /"'o  . f 


t) 


dS^  (38) 


Similar  manipulations  on  Equation  (35)  lead  to  the  following  equation  which 
will  prove  very  useful: 


p(t,tQ,x")  = - ^F(t,tQ,x") 


(39) 


where 


..A 

■"'•'o'-  > •],_ vor;v 


dS 


(40) 


and 


F(t,tQ,x") 


f E„H(x»t  -t)G  (x,x", 

■L  rn?7v- 


t) 


dV 


(41) 


Equation  (32)  shows  that  the  Icernel  of  Equation  (37),  except  for  the  factor 
E^(x,tQ-t)/E^„ (x,tQ-t)  in  the  surface  term,  is  a normalized  kernel. 


11 


The  factor  identically  eqiial  to  unity  if  the  ideal  biasing 

function  is  used.  The  ideal  biasing  function  is  achieved  if  the  "floating  im- 
portance volume"  0^  coincides  with  the  volume  of  the  configuration  (2^  for  all 

X It** 

points  y.  The  factor  is  expected  to  deviate  little  from  unity  is  a™*  ^2 

closely  match  in  the  neighborhood  of  both  x and  x". 

An  integral  representation  similar  to  Equation  (37)  can  be  derived  for  the 

kernel  K (x',x,t  ) defined  by  Equation  (28).  For  that  purpose,  consider  Equation 
^ u 

(36)  with  x"  on  the  common  part  of  euid  E^,  £uid  take  the  normal  derivative  / 

3/3n"  of  that  equation  at  point  x".  Dividing  the  result  by  -Q(x",tQ)  (defined 
by  Equation  10) , one  obtains : 


3 

3n" 


G^(x",x’,tg) 


Q(x",tQ) 


Gq(x',x" 

Q(x",tQ) 


2x(X'Vt) 

. . 3^ 

(x,x*,tQ-t) 

Ex"(x,to-t)  ^_Gq(x,x" 

Ex(x,to-t) 

Q(x-,tQ) 

(42) 


ft) 


dS 


X 


Writing  the  above  equation  for  1>2,  integrating  x*  over  02'  and  taking 
Equation  (10)  into  account,  one  obtains  am  equation  similar  to  Equation  (38): 


L 


e(x",tQ) 


dV 


(43) 


^ Jo  K 


dS 


Equation  (42)  is  an  expression  for  Kg(x",x',tQ)  (Bcjuation  28)  Involving  an 
integral  over  K^(x,x' ,tQ-t)  (Equation  27)  which,  in  turn,  satisfies  the  integral 
Equation  (27) . 
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Finally,  the  equivalent  of  Equations  (39-41)  is: 


1^ 


! 


q(t,tQ,x") 


3^H(t,to,x-) 


(44) 


where 


f 3n"3n 

q(t'to'X")  = / 00^771 


t) 


dS 


(45) 


and 


' ' o'*  ’ -j„  s(x-,t„) 


t) 


dV 


(46) 


Vw2  The  Random  Walk 

V.2.1  Saii5)ling  K^(x",x',tQ) 

The  sampling  problem  is  the  following.  Given  x"  and  t^,  sample  x'  from 
(x",x' ,tQ)/E^„ (x",tQ)  which  satisfies  Equation  (37).  Considering  the  rhs  of 
that  equation,  one  has  to  sanqple  two  terms.  The  first  term  (G^/E)  can  be  san^led 
directly.  To  S2imple  the  second  term,  one  Ccui  first  seui^le  t and  x from  the 
surface  part  of  the  kernel,  and,  given  t and  x,  calculate  the  weight  factor 
E (x,t_-t)/E  „(x,t^-t)  and  sample  G. (x,x' ,t^-t)/E  (x,t^-t).  The  procedure  just 
described  can  be  repeated  for  that  sampling. 

As  shown  in  the  previous  section  (Equation  38),  the  sum  of  the  normalizations 
of  the  density  functions  of  the  first  euid  second  term  is  unity:  only  one  of  the 
two  terms  need  to  be  sampled,  with  probability  equal  to  the  normalization  of  its 
kernel. 

The  normalization  of  the  surface  kernel  is  equal  to: 

where  F is  defined  by  Equation  (39).  With  that  probability,  the  time  variable 
has  to  be  selected  from  the  distribution  p(t,tQ,x")  given  by  Equation  39  or  40. 
Once  t has  been  selected,  the  distribution  of  xcS^  is  proportional  to 
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(47) 


with  remaining  probability  F (tQ,tQ,x'*) , the  volvane  term  of  Equation  (37)  is  to 
be  sampled  for  x*.  The  distribution  of  x'  is  proportional  to 


r(x')  « 


(48) 


Once  x'  is  sampled  from  the  volume  term,  the  random  wal)(  terminates. 

V.2.2  San«)ling  Kg(x",x',tQ) 

The  s^tmpling  problem  is  quite  similar  to  that  described  for  K.  ; given  x" 

¥ 

— 3 

and  t^,  sample  x'  from  y-Tr  G, (x'',x' ,t-)/Q(x",t-)  which  satisfies  Equation  (39). 

0 1 o u 

The  first  term  can  be  sampled  directly.  To  san^le  the  second  term,  one  samples 
t and  X from  the  surface  part  of  the  )cernel,  £md  given  t and  x,  calculates  the 
weight  factor  E^(x,tQ-t)/E^„ (x,tQ-t)  and  samples  K^(x,x' ,tQ-t)  » (x,x' ,tQ-t)/ 

E^(x,tQ-t)  by  the  procedure  outlined  in  the  preceding  subsection. 

Equation  (40)  shows  that  the  sum  of  the  normalization  of  the  two  density 
functions  is  unity:  only  one  of  the  two  terms  need  to  be  sampled  with  appropriate 
probability. 

The  probability  to  sample  the  surface  Icernel  is  equal  to: 


1 - H(tQ,tg,x") 


where  H is  defined  by  Equation  (46).  With  that  prol>abillty  the  time  variable  has 
to  be  saqpled  fzxxa  the  distribution  q(t,tQ,x")  given  by  Equation  (44)  and  (45). 
Once  t has  been  selected,  the  probability  distribution  function  of  xes^^  is  pro- 
portional to 

32 

S(x)  « E^„(x,tQ-t)  GQ(x,x",t)  (49) 


With  remaining  probability  H(t^,tQ,x") , the  volume  term  of  Equation  (39)  has  to 
be  sampled  for  x*.  The  distribution  is  proportional  to 

S(x')  « Go(x',x-,tQ)  (50) 

Once  x'  has  been  sasQ>led  in  this  way,  the  random  wal)c  terminates. 
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VI,  SPECIALIZATION  TO  RECTANGULAR  PARALLELEPIPEDS 


The  derivations  up  to  now  involved  completely  arbitrary  volumes  and 
the  only  restrictions  being  that  the  surface  surrounding  must  be  tangent 
at  a given  point  to  the  surface  of  the  configuration,  and  that  the  volume  £2^ 
must  be  internal  to  both  £2^^  and  £22,  In  practice,  the  choice  of  these  two  volumes 
is  limited  to  such  shapes  for  which  the  Green's  functions  G^  and  G2  are  known  or 
easily  computable,  Carlslaw  and  Jaeger  give  Green's  functions  for  a variety  of 
shapes.  The  efficiency  of  the  Monte  Carlo  technique  would  in5>rove  if  the  shapes 
are  chosen  to  match,  as  closely  as  possible,  the  Ixjundaries  of  the  configuration 
under  investigation.  For  simplicity,  we  limit  the  choice  to  rectangular  parallele- 
pipeds, the  edges  of  £2^^  eind  £22  being  parallel.  This  restriction  permits  an  exact 
solution  in  the  case  of  configurations  with  piece-wise  planar  boundaries,  or 
solutions  to  an  arbitrary  degree  of  accuracy  if  curved  boundaries  are  involved. 

VI. 1 Separation  of  Variables 

We  are  looking  for  the  solutions  G^(x,x',t),  i=0,2,  which  satisfy  Equations 
(4-6).  Let  x'-x  « (Xj^,X2,x^)  in  a coordinate  system  parallel  to  the  cuces  of  the 
RPP.  Dropping  the  subscript  i,  we  are  looking  for  the  solution  G (Xj^,X2,X2,t) 
which  satisfies 


‘-ax,  ax-  ax-  -* 


at 


(51) 


for  i £ a^  , j ■ 1,3  and  t^O  with  the  laoundary  condition 


G(Xj^,X2,x^,t) 


+ aG(x^,X2,X2,t) 

“j 


, - .j  . j 


1,3  (52A) 


+ + 


where  a^  , Oj  are  given  constants: 


Oj  • 0 for  j “ 1,2,  - 0,  « a 


(52B) 
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and  the  boundary  condition 


G(Xj,X2,X2,0)  » 5 (Xj^)  (Xj)  ^ (X3) 


Let  us  assume  that  the  solution  can  be  written  in  the  ionn: 


GCXj^rX^/Xj/t)  “ X^(Xj^,t)  (x^,t) 


Substituting  (54)  into  (51)  we  obtain: 


,x.v  t ^[„5i4.|£].o 

j=l  ^ -> 


A solution  of  (51-53)  is  therefore  (54)  with 


2 

D X^(x,,t)  - X^(x.,t)  = 0 

3Xj  ^ ^ 


+ 3X-’(x.,t) 


* - “j  — 


X-^(Xj,0)  ••  6(Xj) 


j - 1,  2,  3 


The  four-dimensional  problem  (51-54)  has  therefore  l>een  reduced  to  three  two- 
dimensional  problems  (56-58). 

Let  us  now  derive  the  appropriate  expressions  for  the  ftmctions  E and  Q de- 
fined by  Equations  (13)  2md  (11) , respectively. 
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From  the  definition  (13) : 


f 

v’  — 


X ^-a  ^ 


^2~^2 


^2~^2 


, + 
’‘3“®  3 


’‘3-^3 


dXj  G (x^-Xj^»  X2“X2»  Xj'x-jft) 


= E^(x',t)  E^(x^,t)  E^(x^,t) 


(59) 


where 


E^ (x* ,t) 


/-  x'-a . 

/ 


(x’-x,t)dX 


(60) 


Let  us  consider  X'  on  the  surface  of  the  configuration,  the  outer  normal 
pointing  in  the  negative  X^-direction,  From  the  definition  (11)  we  get: 


Q(x',t)  = E^(xj^,t)  E^(x^,t)  Q^(x^,t) 


(61) 


where 


(x’,t) 


-X'-a. 


'x*-a. 


dx 


^ X^(x'-x,t)dX 


(62) 


VI. 2 Sang>ling  the  hPP  Green's  Functions 
VI. 2. ^ For  Purposes  of  Sampling 


Substituting  the  results  of  the  preceding  subsection  into  Equation  (41), 
we  obtain: 


F(t,t„,^")  - F^(t,tQ,x]J)  F^(t,tQ,x|J)  F^(t,tQ,x'^) 


0'  2' 


where 

i /'’‘i"®i  E^(x,t  -t)X^(xV-x, 

F^(t,to)  - -2 i 

Jx^-.a;  ENx^,tQ) 


t) 


dX 


(63) 


(64) 
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Similarly,  liquation  (40)  gives: 


. ■*  123  123  123 

p(t,tQ,x")  = p F p-*  + P p F + F F p^ 


(65) 


where  • 

i ^ ’'‘'“I-''  * I;  ’‘‘‘'■I-*' 

p (t.t^)  - 

E (x^,to) 


(66) 


where 


E^(x,t)  • E^(x,t)  if  X 0 
s 

= 0 if  X = 0 


(66a) 


Both  F^  (Equation  64)  and  p^  (Equation  66)  are  also  functions  of  XV. 
Substituting  63  and  65  into  40,  we  obtain 


P (t,tQ) 


(67) 


Let  us  recall  our  aim:  with  probability  l-F(tQ,t^,x")  semple  the  surface 
term,  i.e.,  first  Scimple  a time  t from  the  distribution  p(t,tjj,x").  We  will 
prove  that  the  following  algorithm  will  produce  such  a time  distribution. 

Perform  the  following  for  i = 1,  2,  3: 

Atten«>t  to  Seunple  t.,  0<t.<t  , from  - |r—  F^(t.,t  ).  This  can  be  done  by 

^ i 

sampling  a random  number  and  attempting  to  solve  F for 

The  atten^t  will  fail  with  probability  F^(tQ,tQ),  in  which  case  set 

When  all  done  (tj^,t2,t2  sampled),  set  t-min  (tj^,t2,t2) . 

Proof:  Given  any  time  T^t^,  the  probability  that  the  algorithm  produces  a 
time  t^T  is  equal  to  the  product  Fj^  (T,tQ) "P^ (T,tQ) 'F^ (T,tQ) , which,  according  to 
Equation  (63)  is  eqxial  to  F(T,tQ,x").  The  probability  that  t<T  is  therefore 
1-F (T,tQ,x"). 

Setting  we  prove  that  the  probability  of  t<tQ  is  equal  to 

l-F(tQ,tQ,x") , as  desired. 

Setting  T»t,  we  can  calculate  the  probability  distribution  function  of 
samples  t delivered  by  the  algorithm 


3t 


[■ 


1-F(t,t 


0' 
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which,  according  to  Equation  (39),  is  equal  to  p(t,tQ,x"),  as  announced.  The 
proof  is  thus  completed. 

Once  t<tQ  has  been  sampled,  the  next  step  is  to  Scunple  x from  the  distribu- 
tion r^(x)  given  by  Equation  (47).  Taking  into  account  separability  of  variables, 
this  can  be  rewritten  as: 


r 

s 


(X) 


= r^(x^)r2(x2)r^X3) 


+ r^(x  )r^(x  )r^(x  ) 
V 1 s 2 V 3 


+ r ' (x  )r  (x  )r  (x  ) 
V 1 V 2 s 3 


(68) 


where 

r^(x)  = E^(x,tQ-t)X^(x^-x,t)/E^(x",tQ)  (69) 

and 

r^(x)  = ENx,t  -t)  ^ X^(xV-x,t)“ 
s s u JX  X 

• ^6(X-X"  + a")  - 6(X-X^  + a^)J/  E^(x",tQ)  (70) 

Expression  (68)  consists  of  the  sum  of  three  terms.  To  sample  (X3,X2,X3) 
from  that  expression,  we  can  sample  a single  one  of  the  three  terms,  with  a 
probability  proportional  to  that  term's  normalization. 

Taking  Equation  (64)  and  (66)  into  account,  we  can  calculate  the  normaliza- 
tion of  the  first  term  of  Equation  (68).  It  is  equal  to: 


p^(t,tQ)  F^(t,tQ)  F^(t,tQ) 


(71) 


The  normalization  of  the  second  and  third  term  of  (68)  can  be  obtained  from  that 
expression  by  a circular  permutation  of  the  indices. 
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Expression  (71)  is  equal  to  the  probability  that  t^^  seunpled  from  p^(t^,tQ) 
is  smaller  than  t^  sampled  from  both  i=2  and  3.  This  property 

can  be  used  to  determine  which  of  the  terms  of  Equation  (68)  is  to  be  s£unpled: 
if  the  time  selection  t=TOin(tj^,t2t2)  produced  t=tj  then  the  term  of  Equation 

(68)  is  to  be  sampled. 

To  sample  the  term,  has  to  be  saunpled  from  r^(Xj^)  (Equation  69), 

for  both  values  of  i/j^  X.  is  to  be  sampled  from  r^ (x  ) (Eqiiation  70);  X is 
^ 3 s j 3 

set  equal  to  ^j~®j  with  a probability  proportional  to 


(72) 


+ 

s ~s  ■ 3 '3'  u " ijx 

Finally,  if  the  volume  term  of  Equation  (38)  is  to  be  sampled,  t is  set  to 


r^  = E^(xV  - a",tQ-t)(+l)  X^  (a;:,t) 


t and  all  x.  are  to  be  sampled  from  r^(x.)  (Equation  69),  for  i*l,2,3. 

0 1 VI 

The  sampling  of  RPP  Green's  functions  has  been  reduced  to  the  sampling  of 
one-dimensional  Green's  functions.  This  sampling  will  be  discussed  in  Section  VIII. 
VI. 2. 2 For  Purposes  of  Sampling  K 


Substituting  the  results  of  Section  VI. 1 into  Equation  (46),  we  obtain  an 
equation  similar  to  Equation  (63) : 


H(t,tQ,x")=  F^(t,tQ)  H^(t,tQ) 

where  F^(t,tQ)  is  defined  by  Equcition  (64)  and 

x;-,;  =hx,t„-t,  ^ X^x-x,t) 


(73) 


3 f ^3”‘ 


(74) 


Q (x'',V 


Similarly,  Equation  (45)  gives  an  equation  similar  to  Equation  (65) 

q(t,tQ,x”)  - p^  F^  + F^  p^  + F^  F^  q^  (75) 
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where  p is  defined  by  Equation  (66)  and 

,2 


q 


E^(x"-a]^,tQ-t)  X^a‘,t)  + E^x-’-aj.tQ-t)  =^X^(a3,t) 

3x 3jc_ 

Q^(x",tQ) 


(76) 


Substituting  (73)  and  (75)  and  (67)  into  (40)  we  obtain  an  equation  similar  to 
Equation  (67) 


q (tjt^)  = 


It  ■'’"'V 


(77) 


As  in  the  preceding  subsection,  the  seunpling  of  time  can  be  achieved  by 
sampling  three  independent  times:  t^  and  t^  from  p^^  and  p^,  respectively,  and 
t^  from  q3(t3,tQ)„  t=tj=min(t^,t2,t3)  is  retained.  If  the  volume  term  is 

to  be  sampled.  Otherwise  the  term  of  the  following  distribution  has  to  be 
sampled; 

rNxi)  shx^)  + r^{x^)  rg(x2)  S^(X3)  + r^(x3)  r^(x2)  slix^)  (78) 

where  r^  and  r^  are  defined  by  Equations  (69)  and  (70)  respectively,  and 


S^(x)  oc  E^(x,tQ-t)  X^(x!^-x,t)/Q^(x^,tQ) 


(79) 


S-'(x)  « 
s 


E3(x,tQ-t)  ^X3(x^-x,t) 

22 

5^x5,  to) 


- r«(x-x5 


3 + a3)-6 (x-x^  + a 


I'] 


(80) 


If  the  term 


is  to  be  sampled,  the  sampling  of  S (x  ) can  be  achieved  by  setting 
+ ® 

X3=X3-a3  with  a probadsility  proportional  to 


(81) 


3i  3 i 2 + 

- E (x^-a*,tQ-t)  (+1)  ^ X^(a",t) 


3X 


(82) 
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Finally,  if  the  volume  term  of  Equation  (42)  is  to  be  sampled,  t is  set  to 

t^,  is  san^led  from  r^(x^)  (Equation  69)  for  i»l,2,  and  x^  is  s^unpled  from 

S^(x,)  (Equation  79) o 
V 3 

The  san^ling  of  RPP  Green's  functions  has  been  reduced  to  the  sampling  of 
one-dimensional  Green’s  functions.  This  sampling  will  be  discussed  in  Section  VIII. 
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VII.  FURTHER  SPECIALIZATION  OF  THE  IMPORTANCE  FUNCTION 


As  defined  by  Equations  (10)  and  (13),  the  importance  functions  Q(x',t) 
and  E^,(x*,t)  are  defined  in  terms  of  a Green's  function  G2(x,x',t),  which,  in 

X 

turn,  is  defined  over  an  RPP  tangent  to  the  outer  surface  of  the  config- 
uration at  the  point  of  closest  approach  to  x'.  We  now  specialize  the  RPP 
to  an  infinite  plane,  as  shown  in  Figure  4. 


Equation  (60)  written  for  j=l,2  gives: 

E^(x|,t)  = E^(x^,t)  * 1 (83) 

Equations  (59)  and  (60)  become: 

E^, (^•,t)  - E^(x^,t)  (84) 

X 

Q(?',t)  « Q^(x^,t)  (85) 

3 3 

In  the  followinr  discussions,  the  superscript  3 of  E and  Q will  be  dropped. 
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The  importance  function  Q is  defined  only  if  the  RPP  Eq  shares  a common 
piece  of  boundary  with  the  surface  of  the  configuration.  The  importance 
function  E does  not  have  to  be  restricted  to  that  case.  We  wish,  however,  to 
inpose  that  restriction,  and  set  E^(x2,t)  ••  1 if  the  RPP  Eq  does  not  abut  the 
surface  This  implies  that  the  importance  biasing  will  be  turned  off  in 

that  case. 


VIII.  THE  ONE-DIMENSIONAL  GREEN'S  FUNCTIONS 


The  sampling  of  RPP  Green's  functions  has  been  reduced  to  the  sampling  of 
one-dimensional  functions  defined  over  the  solutions  of  one-dimensional  differ- 
ential equations  (56).  As  spelled  out  in  Section  VI,  the  time  vari2d)le  has  to 
be  sampled  from  the  p.d.f.  p^(t,tQ)  (Equation  66)  related  to  the  cumulative  distri- 
bution P^(t,tQ)  (Equation  64),  or  from  the  p.d.f.  q^(t,tQ)  (Equation  76)  related 
to  the  cximulative  distribution  H^(t,tQ)  (Equation  74).  The  special  vari6d:>les 
have  to  be  sampled  from  the  p.d.f.  r^(x)  (Equation  69)  or  S^(x)  (Equation  79), 
or  from  the  discrete  distribution  r^  (Equation  72)  or  (Eqxiation  82).  The 

analytical  expression  for  all  the  distribution  functions  involved  will  be  derived 
in  this  section.  Methods  to  Scunple  these  distributions  will  be  discussed  in 
Section  IX, 

1 2 

VIII. 1 The  X , X , and  Related  Functions 

+ 

12  “ 

X and  X satisfy  Equations  56-58  with  » 0.  We  consider  the  symmetric 

case  a^  = -a”  = a/2.  In  that  case,  the  solutions  X^  and  X^  have  been  derived 

previously^.  Summarizing  the  results,  the  solution  can  be  expressed  in  terms  of 

reduced  variaibles 


C = x/a 
T = Dt/a^ 

X^(x,t)dx  = G(C,T)d5 


(86) 

(87) 

(88) 


and  given  either  in  the  eigenfunction  expemsion 

G(C,i)  ■ 2 Z cos  (2n+l)TTC  exp  1- (2n+l) 

n«0  L -I  -I 


(89) 


or  in  the  image  expansion 


' 2/rr  „L 


(90) 
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An  excellent  approximation  to  Equation  (64)  with  E^(x,t)  ■ 1 is  given 
as  follows: 


F(t)  - R(t) 


R(t) 


e ^ du  for  t<t 


1/4 /r 


R(T) 


U' 


exp(-u^T)  - (l/3)exp  (-9Tr^T)  I t>t 


T^T 


T =•  .05  , 

e 


Ola) 

(91b) 

(91c) 

(91d) 


the  relative  error  being  less  than  one  part  in  ten  thousemd. 

In  terms  of  the  reduced  variables  (86)  and  (87),  Equation  (67)  becomes 

P^(t,Tq)  * “ ^ RO)  (91e) 

Equation  (69)  becomes: 

r^(x)dx  « G(5,T)d^  (92) 

Adequate  accuracy  results  from  the  use  of  the  first  two  terms  of  Equation  (89) 
if  t<.05,  amd  of  the  expansion  (90)  if  x^.OS, 

Finally,  expression  (72)  becomes: 

r^  - 1/2  (93) 

s 

meaning  that  x^  is  to  be  set  to  x^  + a/2  with  equal  probability. 

Reference  1 gives  a detailed  description  of  efficient  Honte  Carlo  algorithms 
to  saui^le  the  distribution  (91),  the  differential  distribution  (92)  and  the  dis- 
crete distribution  (93).  Ibese  descriptions  will  not  be  repeated  here. 
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VIII. 2 The  and  Related  Functions  in  the  Case 

The  constant  is  defined  by  Equation  (52b).  It  is  equal  to  the  value 
a(x)  appearing  in  Equation  (6)  and  in  Equation  (3).  The  analysis  up  to  now 
was  completely  general,  assuming  any  non-negative  value  of  a.  Further  analysis 
for  a>0  is  not  completed  yet.  We  restrict,  therefore,  further  discussion  to  the 
case  a^o,  which  applies  to  the  case  of  known  temperature  boundary  conditions  as 
shown  by  Equation  2. 

Let  us  first  calculate  the  importance  functions  E and  Q.  The  solution  of 

^ + 

Equations  56-57  for  j“3,  a^  * - x^,  a^  = “,a  ^ ® 

X^(x,t)  » lim  — — — exp  (-X-X* ) V4Dt)  - exp  (- (x+x* ) V4Dt 

x'-*-0  2/iTDt  *- 

Substituting  that  expression  into  Eqtiation  (60)  one  obtains: 


E(x',t)  = 


s'*  du  = erf  (xV2*'^t) 


x’/2<^ 

Substituting  Equation  (94)  into  (66a)  we  obtain 


E (x',t)  = erf(xV2*^) 
s 

The  subscript  e can  therefore  be  dropped  in  that  case. 

Substituting  (94)  into  (62)  and  letting  x'-H)  one  obtains 


(95) 


(95a) 


Q(0,t) 


1 

/wot 


(96) 


3 3 3- 

VIII. 2.1  The  Functions  p (t,t«,x"),  r (x) , r 

U V s^ 

Let  us  consider  the  symmetric  case  ® j “ ” ®j*  function  X^(x,t)  then 

1 2 

satisfies  the  same  set  of  equations  as  the  function  X and  X discussed  in 
Section  VIII. 1.  In  terms  of  the  reduced  variables  (66)  and  (87)  the  solution 
is  given  by  either  Equation  (89)  or  (90). 


i 

i 
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In  terms  of  reduced  variables,  the  importance  function  E (Equation  (95) 


becomes : 


E(5,t)  =«  erf(c/2v'T) 


where  ^ = 1/2  + t measures  distances  from  the  boiondary.  Let  us  substitute  (97) 
^md  (88)  into  (64) . For  t = Tq  we  obtain 

R(t  ) 

erf  (1/4/r^) 

Ma)cing  the  same  substitutions  in  Equation  (66)  we  obtains 


P (t,Tq) 


erf  (0/2/tq-7)  ■^G(>-  i,T)  + erf  (1/2Aq-t)  GCj^t) 

erf  (1/4/r^) 


The  first  term  veuiishes  as  erf  (0)  • 0„  Ta)cing  (93),  (91a)  and  (67)  into  account, 
the  above  equation  can  be  written  ass 


P (T,T^) 


erf  (1/2/^)  jgR(T) 
erf  (1/4/?^) 


When  R(t)  is  defined  by  Equations  91b,  91c. 

Algorithms  to  seunple  the  time  distribution  (99)  eure  discussed  in  Section  IX. 1. 
Substituting  (97)  and  (88)  into  Equation  (69)  we  obtains 

r^(0  « erf(i;/2/r)  G (?-1/2,t)  (100) 

where  G(C, t)  is  given  by  Equation  (89)  or  by  Equation  (90). 

Algorithms  to  sample  the  special  distribution  (100)  are  discussed  in 
Section  IX. 2. 

Finally,  Equation  (72)  becomes 


0 

s 

The  sanqpling  of  that  discrete  distribution  is  trivial. 
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VIII„2w2  The  Functions  q^Ctjt^)  S^(x), 

The  Green's  function  X^(x^,t)  satisfies  Equations  (56)- (58)  written  for 
j»3,  with  = 0,  a^  = 0 and  a^  = a. 

Let  us  introduce  the  reduced  variables 

C “ *3/2 
T ” Dt/a^ 

X^(x^,t)dx2  “ X(5,T)d^ 


(102) 

(103) 

(104) 


In  terms  of  these  varlabl'  't.  Equations  (56-58)  become: 


X(5,t)  - X(c,t)  «»  0 

(105) 

X(c,t)  • 0 for  ^ * 0,  c • 1 

(106) 

X({;,0)  -=  6(^) 

(107) 

The  solution  of  (105)- (107)  can  be  given  in  the  form  of  an  eigenfunction 
expansion: 

" _ 2 2 
X((;-{;*,t)  “ lim  2 E sin  (mTrc* ) sin  (mire)  ^ ® i 

5*-»0  m*=l 


giving  the  expression 


lim  X(i;-i;',T)  = 


which  converges  rapidly 


2 E mir  sin (miri;)e”™  ^ 
mfl 

for  large  t. 


(108) 
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The  solution  can  also  be  given  in  the  form  of  an  image  expansion 
X(c-i;',T)  ■*  — — |^e;q?(-(C-C')^/4T)  -exp  (-c+c*  )^/4t) 


+ z <-exp(-(c+C'-2n)^/4T)  + exp(-(^-(;'-2n)^/4T) 


+ exp(-(c-c'  + 2n)^/4T)  -exp(-(C+C'  + 2n)V4T) 


giving  the  expression 


lim  ^ 


T X(i;-c',T)  = I 

2T/ivT  ^ 


C exp(-^  /4t) 


+ Z A-(2n-5)  exp{- (2n-c)^/4T)  + (2n  +5)  exp(-(2n  + ^)^/4t)  *| 
n-1  [_  J-* 

which  converges  rapidly  for  small  . 

In  terms  of  the  reduced  variables  (102) > (103),  the  expression  of  the  im- 


portance functions  are: 


E(c,t)  * erf(c/2/r)  "j 

Q{0,t)  - J 


Substituting  (110)  and  either  (108)  or  (109)  into  Equation  (74) , we  obtain 
the  late  and  early  time  expemsion  of  H^(t,tQ).  The  expressions  are  given  in 
Appendix  B.  For  t - t^,  an  excellent  approximation  gives: 


3 "^^‘*^0  ^ 

H (tq/To)  - 1 - 2e  + 2e  for  T<Tg 


(Ilia) 


4/TrT„  e 


for  t>t. 


(111b) 


T - 0.225 
e 


(lllc) 


corresponding  to  a relative  errur  of  less  than  one  part  in  ten  thousand. 
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Let  us  now  rewrite  Equation  (76)  in  terms  of  the  reduced  vari£d3les  (102), 
(103) , taking  expressions  (110)  into  account. 


<3  ^TfTQ)  ® /tttq  erf  (1/2/xq“t)  2 


In  order  to  evaluate  expression  (112) , let  us  start  from  the  late  time  ex- 
pression  of  X(C,t)  (Equation  108)  and  calculate  the  derivative  at  C = 1: 


2 2 2 
. _“,,,n22-nTTT 

—7  x(i;,T)  = - 2 i;  (-1)  n IT  e 


(112a) 


X(c,t)  = ^ 
3? 


where 


R(t)  = -2  r (-1)"  e""  ^ ^ 


If,  instead  of  the  expansion  (112a),  we  start  from  the  early  time  expansion,  we 
obtain: 


[a.^ 


) exp(-l/4T) 


+ I <(1  _ ) exp(-(2n-l)V4T)j 

n=l 


+ (1-  ^ exp(-(2n+l)  V4t) 


r/irT  n*l 


Z (1  - 


(2n-l)  j 


exp(-(2n-l)  /4t) 


an  expression  which  can  be  brought  into  the  form  (113)  with: 


R(t)  - 1 — I exp(-(2n-l)  V4t) 

/ItT 
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Elxpressions  (114)  and  (115)  are  defined  within  addition  of  an  arbitrary 
(and  irrelevant)  constant..  The  constants  have  been  chosen  so  that  expressions 
(114)  and  (115)  are  different  expansions  of  the  seune  function  of  t,  so  that 
R(t)  « 0 at  t = It  also  happens  that  R(0)  * 1.  R(t)  (see  Equation  113)  can 

therefore  be  treated  as  a cumulative  probability  distribution  function. 
Substituting  Equation  (113)  into  (112)  we  obtain: 

q^(T,TQ)  • /itTq  erf  (1/2/Tq-t)  (116) 

where  R(t)  is  given  by  either  Equation  (114)  or  Equation  (115).  An  excellent 
approximation  is  suggested; 


for  T<T 

(117a) 

/ITT 

e 

2 

R(t)  - 2 e"’' 

for  T>T 

e 

(117b) 

T - 0.225 

(117c) 

The  approximation  provides  at  least  four  place  accuracy  for  0_<t<^  . 

Methods  to  seui^le  q^(T»TQ)  of  Equation  (116)  are  discussed  in  Section  IX. 1. 
Writing  Equation  (79)  in  terms  of  the  reduced  varieUsles  (102),  (103),  we 
obtain : 

“ erf  (^2/tq-t)  (118) 

where  X(r,T)  is  given  by  Equation  (108)  or  (109).  We  suggest  that  only  the 

ac 

first  two  terms  of  Equation  (108)  be  leapt  if  t>t^. 

Methods  to  sample  q^Cc)  discussed  in  Section  IX. 2. 
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Finally,  Equation  (82)  reduces  to: 


s 


Sampling  of  that  discrete  distribution  is  trivial^ 


(119) 

(120) 
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IX.  SAMPLING  ALGORITHMS 


IX. 1 Sampling  p^(t,Tq)  and  q^(T,TQ) 


Let  us  define  the  general  function 


p(T,Tg)dT  = N (Tq)  -W^erf  (l/2Ag-T) » ^ R(T)dT 


(121) 


with  a normalization 


)dt  = 1 - N(Tq)R(Tq) 


(122) 


The  function  p(t,Tq)  can  be  made  equal  to  p (T/Tq)  (Equation  99),  by  setting 


N(to)  » 1/erf  (1/4/^) 


W = 1/2 


(123) 

(124) 


and  defining  R(t)  by  Equations  (91b-d) . 

Alternatively,  the  function  p(t,Tq)  can  be  made  equal  to  q^(T,Tg) 
(Equation  116)  by  setting 


N(to)  * H (To,Tg)/R(To) 


(125) 


W =»  /itTq  R(Tq)/H^  (Tg,TQ)  (126) 

and  defining  R(t)  by  Equations  (ll7a-c). 

The  general  sampling  problem  we  discuss  here  is  the  following.  With  proba- 
bility f'(TQ)  set  t>Tq;  else  sample  t from  a renormalized  p(T,Tg). 

Two  different  sampling  algorithms  are  suggested,  each  having  a different 
rzmge  of  efficiency: 
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a.  Small  and  Intermediate  Values  of 

The  following  rejection  technique  is  efficient: 

Step  1 - Sample  t»  0<t<®  from  R(t). 

Ot 

If  t^Tq»  accept  the  S2unple. 

If  t<Tq  do  the  following: 

Step  2 - With  probability  (1-w)  reject  the  sample  and  repeat 
from  Step  1.  With  remaining  probability  w do  the 
following: 

Step  3 - With  probability  erf  (1/2 accept  the  sample. 

With  remaining  probability  reject  the  sample  and 
repeat  from  Step  1. 

The  probability  of  the  algorithm  producing  an  accepted  time  t<Tq  with  dx  at  the 
first  step  is  equal  to 

w erf  (1/2/tq-t)  ^ R(T)dT 

which  is  indeed  proportional  to  the  distribution  (122). 

The  probability  of  rejecting  the  first  Seunple  is  equal  to 


J ° l^l-w-erf  (1/2/tq-t)[]  ^ R(T)dT 


8S  1 


N(xo) 


The  efficiency  of  the  algorithm  is  therefore  equal  to  1/N(tq).  The  proba- 
bility of  sampling  t<Tq  within  dx  is  not  only  proportional,  but  equal  to  the 
distribution  (122). 
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The  efficiency  1/N(Tq)  is  100*  for  'Tq'*’®-  In  both  cases  of  Equations  (125) 
and  (126) , it  is  assymptotically  w/t/irr^  for  large  Tq.  It  becomes  unacceptable 
for  large  values  of  Tq. 

Step  No.  3 involves  a game  of  ch^ulce  with  probability  erf  (1/2Aq-t)  , with 

both  Tf,  and  t given.  The  following  algorithm  can  be  used: 

2 ~ 2 

Sample  a Gaussian  variable  X (p(x)  = — e dx,  0<x<®) . Set  = Tq-1/4x  . 

/tT 

The  probability  that  t>Tj^  is  erf  (1/2/Tq-t)  . 

The  details  of  sampling  in  the  case  of  Equation  (125)  and  (126)  will  be 
given  in  Section  IX. 1.1  and  .2,  respectively, 
b.  Large  and  Intermediate  Values  of 

The  algorithm  just  described  involved  a rejection  technique  based  on  the  in- 
equality erf  (x)  j<l,  which  becomes  assymptotically  an  equality  as  x->«“.  We  now 

2 

propose  to  ta)ce  advantage  of  another  inequality,  erf (x)  ^ — x,  which  becomes 

ZiT 

assymptotically  an  equality  as  x-K).  Having  that  in  mind,  we  rewrite  Equation 
(121)  in  the  form: 


P(t,To) 

= Ni  . fj^  (t)  . f2  (t) 

(127) 

fj^(T)  = 

/^(Tq-t)  erfd/^Ag-T) 

(128) 

fjCl)  = 

f3(T)/M 

(129) 

= 

^ R(t)  t<t„ 

° 

(130a) 

m 

- ^ R(t)  t>t- 

(130b) 

dr  — 0 

r m 

M - 

/ f3(T)dT 

(131) 

Jo 

N ■ 

N(To)/M. 

(132) 
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Similarly,  Equation  (128)  can  be  rewritten  as: 


(133) 


To  sample  Equation  (127),  we  propose  the  following  rejection  technique: 

Step  1 - Sample  t,  0<t<“>,  from  If  accept  the 

sample.  If  t<tQ,  do  the  following: 

Step  2 - With  probability  fj^(T)  accept  the  sample.  Else  reject 
the  sample  and  repeat  from  Step  1. 

The  efficiency  of  the  technique  is  100%  for  > l^ut  deteriorates  for  small 
values  of 

Step  No.  2 involves  a game  of  chance  with  probedDility  fj^{T)  defined  by 
Equation  (128).  In  order  to  construct  an  appropriate  algorithm,  let 

X = 1/4(Tq-t) 

In  terms  of  x,  the  probability  (128)  becomes 


- jji 

Equation  (134)  can  be  expanded  in  Taylor  series: 

2 3 , 

X X X ^ ^ (-x)  ^ 

*^1  “ “ 3 215  " 317  “ “ nl(2n+l) 


(134) 


(135) 


for  x<3  the  absolute  value  of  each  term  of  the  expansion  (135)  is  smaller  than 
that  of  the  preceding  term.  This  property  permits  the  use  of  a particularly 
simple  algorithm: 

Set  n»l,  set  u = a random  number 

Step  1 - Set  u = u - x”/(nl (2n+l)).  Accept  the  sample  if 
u ^0.  If  u<0,  perform  the  next  step. 

Step  2 - Set  u * u + ( (n+1) 1 (2n+2) ) . Reject  the  sample 

if  u<0.  If  it>0,  set  n»n+2  emd  perform  Step  1. 
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If  x^3  (which  is  a rave  event) , we  take  advantage  of  the  semi  convergent 


expansion 


» 1 13 

^ ^ * — * 
(2x) 


1.3.  — (2n^l) 


The  expansion  has  the  property  that,  if  truncated,  the  remainder  is  less  them 
the  cibsolute  value  of  the  first  term  neglected,  and  of  the  same  sign.  This 
property  permits  the  following  algorithms: 


Set  Uq  = a random  number 
Step  1 - Set  u = u 


Step  2 - Set 


Reject  the  sample  if  u<0.  If  u>0,  perform 


the  next  step. 


Step  3 - Set  u = u»x«e  - 1.  Accept  the  sample  if  u^O.  If  u<0, 
set  u = u.2/v^,  set  n = 1 and  perform  the  next  step. 

Step  4 - Set  u = u + (1.3.  (2n-l) )/ (2x)”.  Reject  the  sanple 

if  u^O.  If  n>0,  perform  the  next  step. 

Step  5 - Set  u » u - (1.3.  (2n+l) )/(2x)”^^.  Accept  the  seunple 

if  u^O.  If  u<0,  set  n * n+2  amd  perform  Step  6. 

Step  6 - If  n<X  + 1/2  repeat  Step  4.  If  not,  perform  Step  7. 

Step  7 - Calculate  erf(»^)  by  other  means.  Accept  the  S2unple  if 

Uo<j  J^erf(i^).  Reject  the  sample  otherwise. 

The  test  on  n performed  in  Step  6 corresponds  to  truncation  of  expansion 
(136)  corresponding  to  a minimum  remainder.  Given  x^3,  the  prolsability  of  exe- 
cuting Step  7 (and  therefore  of  having  to  calculate  erf(»^))  is  less  than 
exp  (-3)  •3*5/6'*  < 0.0006. 

The  details  of  san^ling  in  the  case  of  Equations  (123,124)  and  (125,126) 
will  be  given  in  Section  IX. 1.1  and  .2,  respectively. 
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IX.lwl  Details  for  Scunpling 

a„  Small  and  Intermediate  Values  of  Tq 

The  general  method  is  described  in  Section  IX„l.a.  Detailed  algorithms 
for  sampling  ^ R(t)  are  given  in  Section  111^4  of  Reference  1„ 

b.  Large  and  Intermediate  Values  of 

The  general  method  is  described  in  Section  IX.l.b.  We  now  will  work  out 
a detailed  sampling  technique  in  the  particular  case  of  R(t)  as  given  by 
Equation  (91b-c)„  As  defined  in  these  equations,  R(t)  has  different  functional 
forms  for  T>Tg  (x^  is  defined  to  be  equal  to  0.05  by  Ekjuation  91d) , We 

will  assume  tl*at  derive  the  expression  of  f2(x)  (Equation  130)  in  the 

three  cases  0<x<Xg,  ^0*^^* 


a.  The  Case  0<x<x 


Substituting  Equation  (124)  and  (91b)  into  (130a)  we  obtain 

2 


f^  (x)dx  = 


2A( 


° "'1/4/7 


dx 


(137) 


Let  us  perform  the  change  of  varicdsles 


1/4 /t 


where  v * 1/4 /r~ 
e e 


Equation  (137)  becomes 


2 

■ ■ e 


ir^x  Q-x 


(138) 


f^ (v)dv 


dv 


(139) 


which  we  rewrite  in  the  form 


f,{v)dv  = a g h (v)dv 
J e e e 

with 


-2v  V 
0 

h (v)dv  = e 2v  dv 

e e 


(140a) 


(140b) 


(140c) 


(140d) 


V » l/4»^  ; T 1/ (4v  - 4v  )^  (140e) 

e e e 

h^(v)  is  properly  normalized  in  its  range  0<v«»  corresponding  to  and 

0<g  <1. 

— e- 

- B.  The  Case  Tg<T<TQ 

Substituting  Eqxiation  (126)  and  (91c)  into  (130a)  we  obtain 
f (r) i 1 ^ (e-’'^^  - i (141) 

' 2/;rr^)  " ' 

-2p^e”’^^^  (l-3e-®A 
i^o-f 


Let  us  perform  the  change  of  variables  u “ ir/tQ-t  . 
Equation  (141)  becomes 

4 -SiT^t  ^0  “ du 

p(u)du  ■ — (1  - 3e  ) e 

/n 


(142) 
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Performing  the  change  of  variable  u = u^-v,  where  ••  Tr/tQ'Te  » 
Equation  (142)  becomes 


, 4 „ , -STT^t, 

p(v)dv  = — (1  - 3e  ) e dv 

which  we  rewrite  in  the  form 


p(v)dv  = a g h (v)dv 
I i I 


(143a) 


where 


= 4(e 


2 2 
-TT  T -7!  T_ 

- e )/(i/ir  u_) 


(143b) 


g„  = (1  - 3e  ‘)  e 


2 -v(u  -v) 
-OTT  T,  e 


(143c) 


h|j^  (v)dv 


-u  V 
e 

e u dv 
e 


-u 


1 - e 


(143d) 


2,  2 


(143e) 


h^  (v)  is  properly  nonnalized  in  its  range  Ckvcu^  corresponding  to  Tfj>T>T^r  and 
0<g  <1. 


- Y.  The  Case 

0 e 


Substituting  (91c)  into  (130b)  we  obtain 


'3'''  -fif 


(144) 


Ma)cing  the  change  of  variables  v « t-t^,  we  obtain  an  equation  which  we 
write  as 
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f^Cvjdv  = olq  hQ(v)dv 


(145a) 


where 


a 


0 


{145b) 


2 

. 1-3  e"’'  ^ (145c) 

_2 

hQ(v)dv  = e IT  dv  (145d) 

hQ(v)  is  properly  normalized  in  its  range,  and  O^g^l, 

Let  us  now  compute  the  results  of  subsection  a,  B,  and  Y; 
f^(t)  is  given  by  an  expression 

where 

A « and  S^  = for  r*e,*’,  0. 

To  san^le  expression  (146) , one  samples  range  "r"  with  probability  3^ 
(r*e,f,  0).  Given  the  remge  r,  one  san?)le3  h^(v)dv,  performs  the  proper  change 
of  vari^a>le  to  obtain  a time  x,  and  calculates  g^.  With  probability  g^  the 
sample  x is  accepted  as  a valid  san^le  of  ^2^^^  defined  by  Equation  (129).  In 
case  of  rejection,  a new  attenqpt  to  sample  is  made,  starting  from  sampling  the 
range  "r". 

It  happens  that  all  the  distributions  h^(v)dv  are  exponential  which  can 
be  sampled  Isy  standard  methods. 
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IX„1.2  Details  for  Sampling  q (t»Tq) 


a.  Small  and  Intermediate  Values  of  t. 


The  general  method  is  described  in  Section  IX.l.a„  We  will 

detailed  schemes  to  s^unple  — R(T)  with  R(T)  given  by  Equations 

cases  are  to  be  considered  depending  on  t<t  or  t>t  , where  t = 

mm'  m 


- a.  The  Case  0<t<t 


p(T)dT  = - R(T)dT  = - §7  [l e 


dx 


_2_  ^-1/4T 

/vT 


^dx 


Let  u = 1/2/x  ; u = l/2/r~ 
e m 


2 2 

/vj  4d  , -u.,  4,_2,.-u. 

p(u)du  — (u  e )du  = — (2u  -l)e  du 

— du  — 


Let  u = u + w 

® 2 2 
A 2 - -u  - 2u  w-w 

p(w)dw  = — (2u-l  + 4uw+2w)e®  ^ dw 

© © 

✓tt 


Let  V = 2UgW 


p(v)dv  = h^(v)dv 


where 


_1 

4x 


I'm 


"*  (1  + 2x  + 8x 
m m 


g » e 


-T  V 

m 


h (v)  * E p^q. (v)dv 
® v-1 


q^  (v) 

-V 

“ e 

; Pi  - 

(l-2x„)/(l+2x  +8xf) 
m m xn 

qj  (v) 

-V 

D ve 

' P2  • 

4x  /(l+2x„  + 8x^) 
in  Ri  in 

2 

(V) 

- ^ 

2 ® 

; P3  - 

»r^/(l+2T  +8^2) 

™ m m' 

X “ X / (1  + 2x  v) 
m m 


now  work  out 
(117).  Two 
min(x^,XQ). 

(147) 


(148a) 

(148b) 

(148c) 

(148d) 

(148e) 
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Tto  saiqple  defined  by  Equation  (148d),  we  set  v = - log(y)  where  y 
is  the  product  of  1,  2,  or  3 random  numbers  with  probability  Pj^,  p^,  p^,  re- 
spectively. 


- 8.  The  case  t <t 
in 


2 

p(T)dT  “ “ ^ R(T)dT  ” “ 2e**’^  ^d 


(149) 


Let  V ■ T-T 

m 

p(v)dv  = 
where  2 

-It  T 


2 e 

1 


m 


(v)u/ 


h^(v)dv  “ e n^dv 

T ■ T + V 

m 


(150a) 

(150b) 

(150c) 

(150d) 

(150e) 


Recapitulating  the  results,  the  time  distribution  is  written  in  the  form 

A [e^  h^{v)  + e,  g,  hj^(v)] 

where 

A * a + oi„  , 0 “ ®-/A  , 00  ® ®o/A 

To  san^le,  one  first  samples  a range  (e  or  1)  with  probability  g and  g , 
respectively.  Once  the  range  r is  given,  one  scunples  v from  the  appropriate  dis- 
tribution h^(v).  Given  v,  one  accepts  the  sample  with  probal>ility  g^.  In  case  of 
rejection,  the  coo^lete  s2unpling  is  repeated. 

In  practice,  it  was  found  efficient  to  slightly  modify  the  general  technique 

described  in  Section  IX. 1. a.  a is  multiplied  by  w (defined  by  Equation  126b). 

e 

The  game  of  chance  based  on  w described  in  Step  2 of  that  section  can  then  be  by- 
passed. 

£ 

f 

! 
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-A.  Large  and  Intermediate  Values  of  “^q 

The  general  method  is  described  in  Section  IX.l.b.  In  the  case  under  con- 
sideration R(T)  is  defined  by  Equations  (117a-c);  it  has  different  functional 
forms  for  t<t  and  t>t  „ We  will  assume  T Three  cases  are  to  be  considered. 


The  Case  0<t<t 


Substituting  (126b)  and  (117a)  into  (130a)  we  obtain: 


f,(T)dT  = i — -i— 

(Tq-T) 


(151) 


Performing  the  same  operations  as  in  subsection  b,  we  can  rewrite  (Equation  151) 
in  the  form 


f,(T)dT  = a g h (v)dv 
3 e e e 

ere 


(152a) 


a ^ i---  n 

^ 1t/T.(T„-T  ) ^ 


1 + 2t  + 8t 
e e 


] 


-1/4t 


e'  0 e 


(152b) 


h^(v)dv  = q(r)dv 

v.-i- 

2/r  1-/7  /r-i 

e e 

where  q(v)  is  given  by  Equations  (149-150). 


(152c) 


(152d) 

(152e) 
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- 0.  The  Case  t <t<t„ 

e 0 

Substituting  (126b)  and  (117b)  into  (130a),  we  obtain: 

/7t(Tq-t) 

Equation  (153)  is  quite  similar  to  Equation  (141)  of  Section  IX.l.l.b^g. 
Performing  the  scune  operations  as  in  that  section.  Equation  (153)  is  rewritten 
in  the  form: 


f^(t)dt  * g^^  hj  (v)dv 


where 


24’  r 


• e 


(v)dv 


-v(u  -v) 
e 


-u  V 
e , 
e u dv 
e 


-u 


1-e 


2.  2 


- y.  The  Case  t>Tq 


Substituting  (117b)  into  (130b)  we  obtain: 

which  we  write  in  the  form 

P3(t)  - Oq  gp  hQ(v)dv 


(154a) 


(154b) 


(154c) 


(154d) 


(154e) 


(155) 


(156a) 
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where 


■H  Tr 


o„  “ 2e 
0 


^0  “ ^ 


hQ(v)dv  = e TT^dv 


V « t-t. 


{156b) 

(156C) 

(156d) 

(156e) 


Recapitulating  the  results  of  subsections  ct,B,  and  y,  the  distribution 
f^Ct)  is  given  by  expression  (146),  The  sampling  scheme  given  at  the  end  of 
Section  IX. 1.1  applies. 

IX. 2 Sampling  and  S^(i;) 

r^(C)  is  defined  by  Equation  (100).  S^(5)  is  defined  by  Equation  (118). 

Let  us  define  the  general  distribution  function 


p((;)dc  «erf  r(^)dc 


(157) 


p(^)  becomes  equal  to  r2(?)  if  r(C)  becomes  equal  to  G(C-1/2,T)  as  defined  by 
Equation  (89)  or  (90).  p(C)  becomes  equal  to  S^(!;)  if  r{C)  becomes  equal  to 


9x 

T— (^,t)  as  defined  by  Equation  (108)  or  (109), 


Two  different  sampling  algorithms  are  suggested,  each  having  a different 
efficiency  in  different  ranges. 


- a.  Small  and  Intermediate  Values  of  Tq-t 
Step  1,  Sample  {,,  0<C<1,  from  r(5) 

Step  2.  With  probability  erf (^/2/TQ-T)/erf (1/2/Tq-t)  accept  the 

seunple.  With  remaining  probability  reject  the  seimple  and 
start  from  Step  1. 

The  efficiency  of  the  rejection  technique  is  100%  for  but  becomes  poor 

for  large  values  of  Tq-t. 
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I 


The  game  of  chance  in  Step  #2  can  be  implemented  as  follows^  Sample  X, 

2 

0<X<1/2/Tq-t  from  a truncated  Gaussian  (e  ^ dx/erf (1/2/tq-t ) » Accept  the  sample 
C if  ?>  2/tq-t  X. 


- bj  Large  and  Intermediate  Values  of  Tq“T 

In  this  range  of  Tq“T  » we  rewrite  Equation  (157)  in  the  form 
P(^)  « fj^(C) 

where 

fj^(c)  = erf  (i;/2/TQ-T)  /n  (Tq-t)/c  (159) 

auid 


“ ;r(5)dc 


(160) 


As  in  Section  IX.l.b,  we  propose  to  sample  ^ from  f2(c)  and  accept  the 
san^'le  with  probability  fj^(c).  An  efficient  algorithm  to  inclement  the  latter  is 
described  at  the  end  of  Section  IX.l„b„ 

The  efficiency  of  the  rejection  technique  is  100%  as  It  becomes 

poor  for  small  values  of  Iq“I“ 

IX. 2.1  Details  for  Sampling  r^(i;) 

- a.  Small  and  Intermediate  Values  of  Tq~t 

The  general  technique  is  described  in  Section  IX. 2. a.  The  technique  in- 
volves san^ling  r(c)  which  is  equal  to  G(5-l/2,t).  Algoritluns  to  perform  that 
s^unpling  are  descrDsed  in  Section  III. 4 of  reference  1. 

- b.  Large  and  Intermediate  Values  of  Tq-t 

The  general  technique  is  descrilaed  in  Section  IX. 2. b.  It  involves  sampling 
“ ;G(C-l/2,t). 
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As  GK,t)  is  an  even  function  of  K,  the  sampling  of  f 2 ^ performed 

as  follows: 

1.  Sample  , -l/2<^<  1/2,  from  G(C,T) 

2.  Set  ^ = 1/2  + C 

3.  Accept  the  sample  C with  probedsility  Else,  set  C = 1/2  - 
Methods  to  sample  G are  given  in  Section  111^4  of  reference  1„ 

3 

IX„2.2  Details  for  Sampling 
- a„  Small  and  Intermediate  Values  of 

The  general  technique  is  described  in  Section  IXo2.a,  It  involves  sampling 
^ X 

r(r, ) = — (C,t)  as  defined  by  Equation  (108)  or  (109)  „ 


- Uw  Case  of  Early  Times  (t<t  , t =0.225) 

e e 


At  early  times  we  propose  to  use  Equation  (109),  In  terms  of  the  variables 
u = ^/2/r  ; Uq  = 1/2/F  (161) 

Equation  (109)  becomes: 

_u2  “ -(2n  u -u)^ 

r(?)d^  = p(u)du  = 2ue  ^ du  + ^ -2(2n  u -u)e 

— _ *1  ^ 


- (2n  u +u) 

+ 2 (2n  Uq  + u) e du 


or: 


p(u)du  = 2e 


-u 


- - .1 

u - S (u^  - u^)  > 
, n n I 

n=l  J 


du 


(162) 


(163) 


where 

+ -4n  u (n  u +u) 

u = (2n  u,,  + u)  e (164) 


If  u„>l/2  (i.e.,  T<1,  which  is  true  in  our  case  of  small  t) , u - u^  > 0 for 
0 n n — 

all  n.  This  implies 

2 

p(u)du  ^ 2ue  du, 

which  suggests  the  following  rejection  technique: 
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2 

1,  Sample  u,  0<u<Uq,  from  2ue  ^ du.  This  can  be  done  by  scunpling  a 

2 1/2 

random  number  C and  setting  u = [Mod  (-log  (^) , „ 

2„  Let  V = u 

3.  Sample  a random  number  C and  set  w “ C”U 
4„  Set  n = 1 

5.  Calculate  u . If  u >v  jump  to  Step  #8.  Else; 

n n 

6.  If  u^  ^ w the  sample  u is  accepted  (the  sampling  is  completed).  Else; 

7.  Calculate  uj^.  If  u^  < w the  san^le  is  rejected.  Repeat  from  Step  #1, 
Else  jump  to  Step  #9. 

8.  Calculate  u^  and  set  v = v - u + u^.  If  v<w  the  sample  is  rejected. 

n n n 

Repeat  from  Step  #1.  Else; 

9.  Set  n = n+1  and  repeat  from  Step  #5. 

2 

Step  #1  involves  the  sampling  of  2ue  du.  It  remains  to  be  shown  that 

1 “ - + 

Steps  #2-9  correspond  to  an  acceptance  probability  of  1 - — E (u^  - u^) . 

n=l 

Let  us  consider  the  rejection  probeibility  for  each  value  of  n.  If  Steps 
6 and  7 are  executed,  rejection  occurs  if  a random  variable  w,  uniformly  distri- 
buted between  0 cind  u satisfies  u^  < w < u^  (<u) . This  has  probability 

(u~  - u^)/u.  If  Step  No.  8 is  executed,  rejection  occurs  if  the  random  variable 
n n 

w satisfies  v-u  + u^  < w < v(<u).  This  has  also  prob^d3ility  (u  - u^)/u. 

n n n n 

1 CD  _ + 

Summed  over  all  n,  the  rejection  probability  is  indeed  — j;  (u  - u ).  As 

u,nn 

n=l 

the  acceptance  probcdiility  is  eqvial  to  unity  minus  the  rejection  probability,  the 

proof  can  be  considered  as  completed. 

An  instructive  though  lengthy  proof  consists  in  examining  the  acceptance 

probability  for  each  value  of  n.  The  acceptance  probcibility  for  n = 1 is  either 

zero  (if  Step  #8  is  executed),  or  (u-u  )/u  if  Step  #6  is  executed  (as  u < w < u 

n n 

has  that  prob^d>ility) . If  it  is  zero  for  n«l,  it  remains  as  zero  for  the  succeed- 
ing n's,  up  to  axui  including  the  smallest  value  of  n,  for  which 
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V = u =•  z (u  - u ) > u 
0 , n n n„  + 1 

n=l  0 


For  n R + 1,  the  acceptance  probability  is 


u - “n^  + l’ 


(as  + 1 ^ '^0  that  probability). 


For  all  n>r»Q  + 1»  the  acceptance  probability  is 


— (u^  - u ) 

u n-1  n 


(as  u < w < u , has  that  probability) . 
n n-1 

Summed  over  all  n's,  the  acceptance  probability  is 


1 - 1 “ + - 

— (v„  -u  ,)+— r (u,-u) 

u 0 n„  + 1 u n-1  n 

0 n=n  +2 


if  ” - + 

■ = L"  ■ j, 


QED. 


- B.  Case  of  Late  Times  (t>t  ) 

e 

At  late  times,  r ( ^)  is  given  by  Equation  (108).  Keeping  only  the  first 
two  terms,  we  obtain: 


■(C)d5  « 2it  ^ ^ sinirc  + 4e  ^ sin2ir5j  dc 


|^sin(Tr 


5)  + G cos  (ire)  sin  (p^)  Trd^ 


] 


where  c RSe”^^  ^ 


(165) 
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Let  cosCitC)  =•  2x-lv  Equation  (165)  becomes 


r(x)dx  =*  (l-e)dx  + e 2xdx  (166) 

To  Scuttle  (166)  one  can  do  the  following. 

With  probability  (1-e),  x is  set  to  a remdom  number.  With  remaining 
probability,  x is  set  to  the  largest  of  two  random  numbers.  Once  x is  sampled: 

^ =*  cos  ^(2x-l)/Tt  (167) 

- b.  Large  and  Intermediate  Values  of  Tq-t 

The  general  technique  is  described  in  Section  IX. 2. b.  P2(C)  is  proportional 
to  ^r(^),  where  r(5)  has  been  discussed  ^d^ove. 

- a.  Case  of  Early  Times  (t<t  ) 

e 

The  early  time  behavior  of  r(^)  has  been  discussed  in  section  IX. 2. 2. a. a. 
Performing  the  change  of  varied>le  u * 5/2/r  we  obtain: 

r(5)d5  <x  u p(u}du 

where  p(u)  is  given  by  Equation  (163). 

The  technique  described  to  s£uiiple  p(u)  can  be  easily  modified  to  sample 
u p(u)du.  Only  the  first  step  of  the  rejection  techniques  needs  to  be  modified. 

_u2  2 -u^ 

Instead  of  sampling  ue  du,  one  samples  sue  du,  which  C2in  be  achieved  by 
setting 

u « J-log(Cj^)  sin^(irt2) 
emd  accepting  the  saii9>le  of 

The  remainder  of  the  rejection  technique  applies  without  modification. 

- 8.  Case  of  Late  Times  (t>t  ) 

e 

At  late  times,  we  propose  a simple  rejection  technique: 

1.  Sample  r(c)  as  discussed  in  Section  IX. 2. 2. a. 8,  0<c<l. 

2.  Accept  with  probability  t;.  Else  reject  emd  repeat  Step  1. 
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APPENDIX  A 


A Crucial  Inequality 


Let  be  a volume  surrounded  by  surface 
Let  x"  be  a point  on  (See 

Figure  6 ) . Let  fi  ^ be  a volume  surrounded 
by  surface  such  that  and  E^  are 
tangent  at  x",  with  the  same  outer  normal. 
Let  Sj  be  the  surface  coirmon  to  E^^  and  E^ 
(S2  can  degenerate  to  the  single  point  x") 
And  let  = Ej^  for  i=l  eind  i*3. 


Let  and  be  Green's  functions  satisfying  Equations  (4,5,6)  for  i=l  and  3, 
respectively. 

We  are  going  to  prove  the  following  inequality: 


Gj^  (x",x,t)dV^ 


G3(x",x,t)dV^ 


< 00 

for  0<t<tQ<a> 


(Al) 


Let  us  introduce  a surface 
surrounding  a volume  0^,  S 
pletely  internal  to  both  E^^ 
also  Introduce  a surface  E^ 
surrounding  a volume  S 

external  to  both  Ej^  and  Z^. 


h,0  " "2 


1,0 

and  E 


is  com- 

Let  us 


1,2 


^1,2  ^2 

is  completely 


FIGURE  6 
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Equation  (42)  of  Section  IV  has  been  derived  without  assuming  (Al).  Let 


¥ 


us  consider  that  equation  for  i“l  emd  3,  and  integrate  both  sides  over  the 

volume  0..  We  obtain: 

1 


Q(x",tg)  Jjj  Q(x",tQ)  x* 


E. 

1 


E 

X 


(X,tQ-t) 

n(X,to-t) 


Ex-'<*'^o"^^  3n3n"  GQ(x,x",t) 
Q(x",tQ) 


dS 


(A2) 


where 


2^(x,t)  G^(x,x’,t)dV^, 


(A3) 


(A4) 


As  the  solutions  of  the  set  of  equations  (4-6)  2u:e  positive,  we  have 


G^(x,x',tQ)  ^0  for  Xjx'eS)^. 

This  implies 
E^(x,t)  >0. 

From  the  )3oundary  condition  (6)  we  derive 

^G^(x,x',t)  («.«•.  V i" 

which  issues 

Q^(x,t)  ^0 

FurtheiTBore,  we  can  show  that  for  t««> 

Qj^(x,t)  >0  for  t«». 


S4 


(A3a) 


(A4a) 


(A4b) 


i 

* 


According  to  Equation  (10)  eind  (A3)  we  also  have 


Q(x",tQ)  - Qjtx-.tp) 


(A5) 


and,  according  to  Equation  (13)  and  (a4) : 


EjjH^Xft)  — E2(x,t) 


From  (A2)  we  can  derive  the  upper  bound 


2i<^"'V 

Q(x",tQ)  - 


A + B‘V(Ej^/E2) 


and  the  lower  bound 


(A6) 


(A7) 


®3^^  ,t  ) 


where 


dV  , 

x' 


r r 3" 

f^o  . f 

Jo  Jh,o 


t) 


dS 


(A8) 


(A9) 


(AlO) 


Equation  (43)  giving 
A + B « 1 


(All) 


U(E^/E^)  and  L(E^/E^)  are  respectively  the  upper  bound  and  the  lower  bound 
of  the  ratio 

(x,t)/E^  (x,t)  for  X on  Q and  0<_t<_tQ. 
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We  now  proceed  to  show  that  both  A and  B are  non- negative,  a result  necess- 
ary for  the  assertion  (A7), 

Using  the  symmetry  property  of  Green's  functions 

G^(x,x',t)  « G^(x',x,t) 


we  can  rewrite  A in  the  form 


Q2(x",V 


QoU-.tfl) 

Q2(x-,V 


T2diing  (A4a)  into  account,  we  obtain 


or,  taking  (A4b)  into  account: 


A>0  for  tQ<“’ 


Now  consider  the  expression  (AlO)  of  B.  Both  E^„(x,tQ-t)  and  Q(x",tQ) 
are  non-negative  (Equ^^tion  A3a  and  A4a). 

The  remaining  term  will  be  shown  to  be  also  non-negative. 

Indeed,  it  can  be  rewritten  as 


Jo(x,x-,t)  - ^ [^GQ(x,x-,t)] 


or,  using  the  boundary  condition  (6) 


Using  the  property  of  r’lnmetry 


Involving  Equation  (6)  again,  we  obtain: 

All  terms  of  the  integrand  of  B (Equation  10)  being  non-negative,  we  have 


(AlOa) 
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In  order  to  determine  the  bounds  of  E^/E^,  we  turn  to  consider  Equation  (36) 
which  we  rewrite  in  the  form 


G.  (x",x',tQ)  = G.(x',x",t  ) 


t)dS 


(A12) 


where  ft  . CS2 . . 

1 D 


(A13) 


Integrating  Equation  (Al2)  over  x'efi.,  we  obtain 

3 


E.(x",to)  = 


t)dS 


(A14) 


which  leads  to  the  inequality 


E. (x,t)  < Ej (x,t) 


(A15) 


provided  (A13)  is  true. 


Let  i-1  and  j=2.  As  (AlS)  gives 


E^(x,t)  ^ E2(x,t) 


or 


W(Ej^/E2)  - 1. 


(A16) 


Substituting  the  above  expressions  into  (A7)  and  ta)cing  (All)  into  account, 
we  obtain 


Ql(x",tQ)/Q(x",to)  11 


(A17) 
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Let  us  write  the  inequality  (A4a)  for  i^l; 


E3(x,t)  >0 


giving 

LCE^/E^)  = 0 


(A18) 


Substituting  that  equation  into  (A8)  we  obtain 


Q^{x”,t^)/Q{x\t^)  > A 

Finally,  Equation  (A17)  and  (a19)  imply 

Ql(x",to)  ^ 

Q3(x"'to)  - 


(A19) 


(A20) 


As  A>0  for  tQ<“’  (see  Equation  A9b) , Equation  (A20)  implies: 


Ql(x",V 

Q3(x",to) 


<® 


(A21) 


Substituting  Equation  (A3)  for  i=l  and  3 into  (A21)  we  obtain  (Al)  Q.E.D. 
The  current  proof  required  only  tangency  of  surfaces  and  The 

introduction  of  a very  special  surface  such  that 

construction  of  the  proof  only.  In  the  rest  of  the  text,  the  restrictions  on 
I 2 are  only  those  imposed  here  on  ^3. 
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APPENDIX  B 


Expressions  for  H (tjt^) 


H (t,tQ)  is  defined  by  Equation  (74)  of  Section  VI„2.2.  Starting  from  the 
early  time  expansion  Equation  (109) , we  derive* 


H^(t,to)  = ^ S2  -~erf(u} 


4t 


m=l 


- e 


"^0  r 1 

^erf  (Uj^+U2)  -erf  (Uj^-U2)J 


where 


J—0 


= i / ^0  _ 1 


starting  from  the  eigenvalue  expansion  Equation  (108): 


3 “>  r „ 2 2.  -m^Jr^t„  -| 

H (t,t  ) = 2/^  Z erf  (u)  (-dV”  ^ -e  ° R 

in=l 


where 


u •=  1/2  '^t^-t 


and 


or 


• j [erf(u  + i 2J)  t erf(u-i  21)] 

_ 2 

‘'m  ■ ^ hr 

[]-(-l)"  Cosh  2^] 


y,  2 „ -n  /4 

I f 

^ n=l  u^+4n^ 


(Bl) 


(B2) 


(B3) 


(B4) 


The  derivations  are  given  in  MAGI's  internal  memorandum  P-7133,  Sept.  13,  1976, 
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